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In a generalized framework, where multi-state and inter-state linkages are allowed, we derive a 
sufficient condition for the stability of synchronization in a network of chaotic attractors. This condi- 
tion explicitly relates the network structure and the local and coupling dynamics to synchronization 
stability. For large Erdos-Renyi networks, the obtained condition is translated into a lower bound on 
the probability of stability of synchrony. Our results show that the probability of stability quickly 
increases as the randomness crosses a threshold which for large networks is inversely proportional 
to the network size. 
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INTRODUCTION 



The synchronization of complex dynamical systems is one of the most intriguing problems and has roots in physics, 
biology, and engineering pQ. The problem of synchronization in a network of interacting nodes was first brought 
to attention by Wiener [21] and later pursued by Winfree, who modeled biological oscillators as phase oscillators, 
neglecting the amplitude [2"][T5]. 

Following Winfree's pioneering work, there has been considerable effort to study the structural and dynamical 
effects of the network and its nodes on the state of synchrony. One important problem encountered in these studies 
is the possibility of separating the global (network) and local (node) characteristics to allow a general assessment 
on stability of the synchronous state. In the seminal work by Pecora and Carroll [TO], this issue was addressed by 
introducing the concept of master stability equation. This equation provides a condition on stability of the network, 
known as master stability condition (MSC). The MSC also leads to useful bounds on structural properties of the 
network such as eigenvalues of the Laplacian matrix of the network [2D] [ID]. Following this important result, most 
of the focus has been directed toward interpreting the bounds on eigenvalues of the Laplacian matrix of the network 
implied by the MSC as bounds on the degrees of nodes in the network. In [5] and [T3], the effect of heterogeneity 
and smallness of the network on synchronization stability have been investigated and it has been shown that if the 
network is more homogeneous the stability of synchrony is easier to achieve. Though most efforts have focused on 
the unweighted and undirected graphs, [7] and |16j have investigated the master stability equation in directed and 
weighted structures, respectively. 

Since the synchronization manifold can be considered as a fixed point of a reduced (induced) space, the vast majority 
of the existing literature on the synchronization in the networks considers the largest traversal Lyapunov exponent 
as a measure of stability of the synchronization. Of course, as noted in |10j . the negativity of Lyapunov exponents 
is not a necessary nor a sufficient condition on the stability of manifold itself and it does not stop the manifold from 
bubbling and bursting [18] [5]. 

In this paper, we use an alternative master stability equation derived from Lyapunov direct method to obtain a 
condition on the stability of the whole network, based on the eigenvalues of symmetric part of local and coupling 
dynamics. Since this condition encompasses the Lyapunov spectrum, it is a sufficient condition on the stability 
[18j . Furthermore, we generalize the conventional setup where the linkage matrices are diagonal (often with binary 
components) to the case of an arbitrary linkage matrix, allowing multi-state and cross-state linkages, possibly with 
different strengths, which is now receiving more attention [I] [5]. 

We then use the derived MSC to calculate a lower bound on the probability of stability for large Erdos-Renyi 
networks. We relate the condition of stability to dynamical characteristics of individual nodes and their coupling and 
structural properties of the Erdos-Renyi networks, namely, network size and randomness parameter. 



where Xj € E™, i = l,...,N, denotes the state vector of node i, and JF(-) and %{■) denote the node and coupling 
dynamics, respectively. L = [l^j] is the Laplacian matrix of the network. Since L has zero row-sum, this network 
has a synchronization state, xo, which is the solution of local state equation, Xo = .F(xo), and xo(0) = Xoo- This 
equation also defines the synchronization manifold. To maintain the synchrony throughout the network, all x^ should 
converge to a synchronous state, x . This means that all the modes traverse to the synchronization manifold should 
be damped out [TO] . Denote the deviation of the state of each node from the synchronous state by = x^ — xo . If 
are small, can be linearized around xo as 



SYSTEM MODEL 



Consider a network of N identical nodes with identical coupling dynamics 



N 




(1) 



JV 




(2) 



where F and H are Jacobian matrices of J-(-) and W(-) around Xo, respectively. Note that in general, Xo is time 
dependent and so are the Jacobians. Here, for brevity of the presentation, we have dropped the explicit dependence 
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of F and H on Xo (and therefore on time) from our notation. Combining ([2]) for all i, we obtain the unified dynamical 
equation of the whole network as 

e = ( Ijy ® F — L (8 H )e, (3) 
F 

where Ijy is the identity matrix of size N, <S> denotes the Kronecker product, and e = [ej e% ■ ■ ■ eJ { } T is the error 
vector with respect to Iat <8> x.o, where ljv = [1 1 • • • 1]jvx1- Superscript T denotes the transpose operator. 

In the rest of the paper, we assume that the network is undirected, and therefore its Laplacian matrix is symmetric, 
i.e., hj — Iji for all i,j. 



NETWORK STABILITY CONDITION 

The network synchronization is exponentially stable around Xo, if F is exponentially stable. Due to zero row sum 
property of L, its smallest eigenvalue, hn, is zero [BJ. The mode for fi^ = in ([3]) represents the variation parallel to 
manifold and hence should be omitted in the study of stability for traversal exponents |10j . Furthermore, multiplicity 
of zero in the set of network eigenvalues is one, if and only if the the graph is connected [BJ . The remaining modes in 
([3]) are traverse to the synchronization manifold and they decay exponentially if and only if the corresponding modes 
of the system governed by F are exponentially stable. 

Now, considering the system ([3]), the norm of e can be bounded by ([3] Thms. 1 and 2) 



i /"' Ai( F(r)+F (t) )+e dr 

- T 

for all t > to and some e > 0. Here Ai(.) denotes the largest eigenvalue of the argument. We note that if F+F +el -< 
for all t > to, the exponent approaches negative infinity as t approaches infinity and the system is stable. Thus, a 
sufficient condition on stability is 

3e > s.t. F(<) + F T {t) + el -< 0, Vt > t , 

here -< denotes negative definiteness of a matrix. This condition guaranties damping of traverse modes and can be 
used as stability criteria in the study of dynamical systems [TU] , 

Since F can be written in block Jordan form, the equivalent sufficient condition on asymptotic stability for ^ 
yields 

F + F T -^(H + H T ) +el^0, Vt > to, (4) 

for 1 < i < N — 1. Hence, F is exponentially stable if and only if F — /XjH is stable for alH — 1, 2, ■ • • , N — 1, where 
fii is the ith largest eigenvalue of L (see appendix for propf). The condition in Q requires the largest eigenvalue 
of symmetric part of F — /ijH to be strictly negative, and not approaching zero from below. Here we have used the 
assumption that L is symmetric. Thus, \Xi are real and non-negative [? ]. In the case of complex F and H, Hermitian 
operator replaces the transpose operator. 

Now we can use Weyl's inequalities to further explore ([4]). The Weyl inequalities provide upper and lower bounds 
on the eigenvalues of sum of Hermitian matrices. Assume that B and C are n x n Hermitian matrices, and 1 < 
j,k,j + k — n< n, then Weyl's inequalities state that ([11] Thm. 1.3) 

A j+fe _i(B + C) < A,(B) + A fc (C), (5) 

and 

X 3+k - n (B + C) > A,(B) + Afc(C), (6) 

where Afc(.) denotes the fcth largest eigenvalue of the argument. Utilizing ^ with j + k — 1 = 1 and ^ with 
j + k — n = 1 yields 

A 1 (B + C)<Ai(B) + Ai(C), (7) 
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and 



Ai(B + C) > A i (B) + A n _ i+1 (C), 



for 1 < j < n. Now, since F + F 1 and -/i, t (H + H J ) are Hermitian, we can use (|7|) and (|8| to bound Ai(F + F — 
/iiH — /iiH T ) + e from both sides: 



(8) 



Ai(F + F T - Mi H - Ml H T ) + e < Ai(F + F T ) + Ai(-^H - Mi H T ) + e 

= Ai (F + F T ) - MiA n (H + H T ) + e, 



(9) 



and 



Ai(F + F T — W H - /i. t H T ) + e > A,(F + F T ) + A„_ J+ i(-^H - ^H 2 



A,(F 



Ml A,(H + H T )+e, 



(10) 



for 1 < j < n. 

To obtain a sufficient condition on stability we now force the upper bound in |9]), to be negative 

Ai(F + F T ) - mm{/XiA„(H + H T )} + e < 0. 



(11) 



for i = 1, • • • , AT — 1 and some e > 0. Since A„(H + H T ) maybe positive or negative, the second term in ( 11 ) reduces 
to min{/iAr_iA„(H + H T ),/iiA„(H + H T )}. Thus, a sufficient condition on stability is 



[i N -i > M f if A„(H + H T ) > 
Aii < /x+ if A„(H + H T ) < 



(12) 



for some e > 0, where ^ = Ai(F + F T ) + e/A„(H + H T ). 

We can also derive a necessary condition for Q by forcing the lower bound in ( 10 ) to be negative, or 

> max |Aj(F + F T ) - ^Aj(H + H T )| + e. 



(13) 



Let k be the index of smallest positive eigenvalue of H + H . Then (131 reduces to 

Mi < Mi! 



(14) 



where 



t A 
Mjv-i — 



A, F + F J 



A, ( H + H 



(h + h t ) 



and 



,t A 



Xj ( F + F J 



J>fc aJh + h 



In the case that Xj(R + H T ) = 0, if Aj(F + F T ) > then the condition Q is not satisfied, and if A 3 -(F + F T ) < 0, 
the corresponding condition can be eliminated. 

We note that, (14) is a necessary condition on Q, which itself is a sufficient condition on stability. Therefore (14) 
does not reveal anything about the stability of the network. However, since (12) and (14) sandwich Q, (14) provides 
some information regarding how close ^ and (fl2| ) are. 

Having developed the bounds on condition ([4j), namely ( 12 ) and (14), we now proceed to relate them to the degree 
properties of the network. To do this, we employ following inequalities known for symmetric Laplacian matrices [5] 



< 



N 



' dm i r 
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and 



Mi < 



N 



where d max and d m i n denote the maximum and minimum node degrees, respectively. Using these, and (12), the 
stability condition can be also expressed as 



< dt 



if A„(H- 
if A n (H - 



H T ) > 
H T ) < 



where 



and 



Similarly, necessary conditions for Q become 



where 



and, 



Jt it 1 t 

a min — M/V- 



dt A*-l..t 

max 



N 



Mi- 



(15) 



tax ^ ^max; 



4 — 



N - 1 



iV 



iV ■ 



N 



: m1, 



From (12) one can draw the conclusion that in the network with low algebraic connectivity, (J-n-i, synchronization 
is difficult to achieve. This behavior is caused by the fact that the coupling is not strong enough to push/pull the 
oscillators to synchronous state. And from ( 15 ) we can see the other case of non-synchronization behavior occurs 
when (some) nodes have too many connections (condition on /ii). This phenomenon is known as synchronization 
quenching, where the coupling is so strong that eliminates the self-drive of (some of) the oscillators and consequently, 
the network cannot achieve synchrony [9]. 



PROBABILITY OF STABILITY FOR ERDOS-RENYI NETWORKS 



In the following, we investigate probability of stability of Erdos-Renyi networks |17) . For large Erdos-Renyi networks 



with randomness parameter p, we can use ( 12 ) and ( 14 ) to calculate the lower and upper bounds on the probability of 
Q being satisfied. We recall that (12) provides a lower bound on the probability of stability, whereas (14) describes 
the closeness of this lower bound and the probability of Q. 

Since eigenvalues of any large randomly generated symmetric matrix follows the Wigner's semi-circular distribution 
[14j . we can approximate the distribution of eigenvalues of Laplacian for an Erdos-Renyi network. By definition 
L = diag(d) — A, where A is the adjacency matrix of the network, and d is the degree sequence of the nodes. Also 
in large Erdos-Renyi network, we can approximate diag(d) by Npl, the average degree of the network [T7]. Hence, 
L R3 Npl — A. Thus the eigenvalues of L have (approximately) the following distribution [T4"] 







Vl-X 2 \X\ < 1 
elsewhere 



where 



X 



x — Np 



2y/Np(l-p) 



G 



The order statistics /ijv-i an d Hi have densities 

p^_ 1 (x) = (N-l)[l-P,(x)} N -%(x), 

and 

p lll (x) = (N-l)[P IM (x)] N - 2 PlJt (x), 

where 



J(N-l)p-2y/N P (l-p) 

= 1 - i cos" 1 X + -Xy/l -X 2 



Now, we can evaluate the probability of occurrence of each conditions given in (12) and (14). The lower bound on 



the probability of stability attained from Wigner's approximation for (|12|) is 

P W,LB 



= \l-[l-Pjrf)] N - 1 A„(H + H T )>0 
\ [P,(^)] A„(H + H T )<0 



(16) 



Similarly an upper bound on the probability of Q being satisfied can be derived by applying Wigner's approximation 
to pi 



AV,UB = 2 (l- [1- PM-i)]"- 1 + <~ Vi) [PM) PM-X' 1 



(17) 



Since these probabilities have relatively sharp roll-offs as a function of N (see numerical results), we can use 
randomness values, p\^ and p\j, which yield Pyy lb(-Pl) = 1/2 and Pyy tjbCptj) = V^i to study the synchrony trends 
as network parameters change. From (16) and (|17|, we have 



1 



PL 



and 



PU 



iV + 4 



1 

N + 4 



Ai(F + F T ) + 6 
A„(H + H T ) 



max 



X^F + F 1 



j>k \j (H + H ) 

Thus, the value of randomness, p, required to have stable synchronous state decreases as 0(1/N). 



NUMERICAL RESULTS 

For a numerical example, we consider a network of Rossler oscillators (with parameters a\ — 0.165, a% = 0.2, and 
03 = 10) [T2] coupled through all of their states. This set of parameters results in a coherent chaotic oscillation, since 
a\ < 0.21 [5]. Therefore, the Jacobian of the oscillator can be computed as 

-1 -1 

1 0.165 

X3 Xi— 10 



We also assume 



H 



V6 + V3 + V2 



V2 V2 V2 

o Vs Vs 

y/6 



where c = trace(H), is the coupling strength. 
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FIG. 1. Probability of stability as a function of p. 



With this setting we calculate our results on the probability of stability of the network for different values of N, 
p, and c. To this end, we have considered N trajectories starting from initial point [sj Sj 0) T where Sj's are selected 
uniformly from interval [0.95, 1.05] and all the corresponding eigenvalues are calculated over average of 20 cycles of 
initiated trajectories. 

Fig. [T] shows the probability of the stability of the network as a function of network randomness, p, for different 
network size, N, and with coupling strength c = 1. As it can be seen, probabilities of ( 12 ) and ( 14 1 are close to that of 
Q. Moreover, we observe that the approximated probabilities provided by the Wigner's distribution of eigenvalues 
of the network are also reasonably close. 

As Fig. [l] shows for positive definite H + H T in a large network if the average degree, pN, is above some threshold, 
PlN w Ai(F + F T )/A„(H + H T ), (in this example approximately 50) the network becomes stable. Note that in 
this particular numerical example, due to positive definiteness of H + H T , only the transition from asynchrony to 
synchrony is observed. 

The behavior of the network in the sense of its stability versus network size for several values of p is shown in Fig. 
[2j Once again we observe that the probability of stability suddenly increases as pN crosses the threshold above. 

Fig. [3] shows that the stability of the network grows as the network size and coupling factor increase. Of course this 
is due to our choice of coupling dynamics which is positive definite and by increasing coupling strength, it provides 
stronger negative feedback to stabilize the network. 



CONCLUSION 



In conclusion, considering an alternative master stability condition, we have derived a sufficient condition of stability 
which is a function of the eigenvalues of network structure and symmetric parts of linearized local and coupling 
dynamics. Our condition relates the largest eigenvalues of the symmetric coupling and the symmetric local dynamics 
to stability conditions of the networks. For Erdos-Renyi network we have calculated a lower bound on the probability 
of stability. Then we have proceeded to calculated associated threshold value of randomness where the system starts 
to become stable as p increases beyond this threshold. The reason for this phenomenon is that below the certain 
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FIG. 2. Probability of stability as a function of N. 
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FIG. 3. Probability of stability as a function of c. 
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threshold, all or some of the nodes cannot achieve sufficient information exchange. As a result, those nodes cannot 
synchronize themselves with the rest of the network. 

Appendix: stability of e = Fe 

Since Lis real and symmetric, it is unitarily and orthogonally diagonalizable (Ch. 5.4. [15]). That is, L = QAQ T , 
where Q is unitary and A is diagonal [T5]. Define e = (Q T ® I)e. Since e are mapped to e by a non-singular linear 
transformation (Q T is unitary), their stabilities are equivalent. We have 

e = (Q T ®I)F(Q(g)I)e, (18) 

using the properties of the Kronecker product we have 

e = (Q T ® I) (I ® F)(Q ® I)e + (Q T ® I)(QJQ T ® H)(Q ® I)e 
= (Q T ® F)(Q <g> I)e + (JQ T ® H)(Q ® I)e 
= (I ® F + A ® H) e. 

The matrix I ® F is block diagonal. Thus the resultant diagonal block matrix has F — /^H as its diagonal blocks. In 
other words, stability of F is equivalent to stability of F — /i.;H for i = 1, • • • , N. 
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